On numerical averaging of the conductivity coefficient using 

two-scale extensions 



Vsevolod Laptev* 



Abstract 

' In this article we compare solutions to elliptic problems having rapidly oscillated conductivity 

, (permeability, etc) coefficient with solutions to corresponding homogenized problems obtained 

^SJ ' from two-scale extensions of the initial coefficient. The comparison is done numerically on 

several one and two dimensional test problems with randomly generated coefficients for different 
' intensities of oscillation. The dependency of the approximation error on the size of averaging is 

, investigated. 
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^ ■ 1 Introduction 

-4— » 

We consider a second order elliptic equation with a rapidly oscillated coefficient aA/(-): 

-V-{aMix)Vu)=f in 17, u\9n = g. (1) 

' The equation appears to describe sucli problems as the stationary heat transfer in composite mate- 

rs ■ rials, the flow in non-homogeneous porous media as well as in many others. For the periodic coeffi- 

! cient, the averaging procedure is well known and is called the periodic homogenization PP , [13] , [IZ] • 

The general non-periodic case is very important for practical applications (e.g. in geoscience, 
petroleum engineering), and a vast literature exist discussing and comparing algorithms intended 
. for averaging the permeability coefficient (see the reviews [5].|llj.[T5].|16j). Some algorithms are 

based on the idea that the effective (averaged, upscaled, equivalent grid block) permeability field 
in the whole domain can be determined by solving the flow problem locally [3]. They vary in the 
^ , choices of the local subdomain, the boundary conditions, and the ways to extract the effective 

^ I permeability coefficient from the solution of the local problem. These algorithms usually perform 

well, and intuitively there should be arguments to justify their usage. 
An effective coefficient A{-) in the averaged problem 

-V-{A{x)VU) = f inn, U\9n=9- (2) 

is generally different from aMi'), and the solution U is different from u. Therefore there is a 
trouble with perfect justification: the difference between these two solutions in some cases can be 
unacceptable. Homogenization is known as a rigorous way to justify an averaging process. This 
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is so because homogenization deals with sequences, not with single problems. And if the sequence 
of problems converges in some sense to a limit problem then, whatever strict requirements we 
have, there is always a set of problems from the sequence for which this limit can be considered 
as an averaged problem. Although in practice we usually need to upscale a single problem like 
(dl), not the whole sequence. Nevertheless, if our initial problem ([TJ belongs to the sequence in 
the homogenization process then the limit problem may be a reasonable candidate for upscaled 
initial problem, even if we cannot improve the approximation. We only need that the sequence is 
homogeneous in the sense that all its members, including our initial problem, have something in 
common (it is important to avoid situations when a convergent sequence contains an element which 
has nothing to do with the rest of the sequence). 

One way to do so is to use the sequence from locally periodic homogenization (see e.g.[r, p. 71]) 

-V-(^a(^x,^^Vue^=f inn, Ue\dn = 9, (3) 

where the function a{x,y) is a two-scale extension of the initial coefficient aM(')' 

Definition 1.1 (from [S]). Let us say that a function a{x,y), {x,y) G $7 x W^, 1-periodic in the 
variable y, is a two-scale extension for ani{x) if there exists a positive number e such that 

a = aM{x), Vx € 0. (4) 

Having a two-scale extension, we can choose a strictly positive sequence {e„} — > 0, containing e, 
and consider ([3]) with e from {£n} as a sequence in the scope of locally periodic homogenization. The 
expressions for A(x) and corrections of U can be found in the literature devoted to homogenization. 
All the members in the sequence ([3]) have in common the function a{x, y); and at e = e we recover 
the initial problem ([1]). In this sense a two-scale extension establishes a connection between ([2]) 
calculated from the homogenization algorithm, and the initial problem ([T]). As it was already 
mentioned, we cannot claim that ([2]) with such A{-) is the averaged problem for ([1]). Moreover, there 
are (infinitely) many two-scale extensions leading to different A{-) for the same aM(')- Nevertheless, 
we expect that among them there could be classes of extensions appropriate for averaging. Therefore 
it is interesting to test numerically the two-scale extensions from [S] on several model problems with 
non-periodic coefficients. In each test we calculate both the solution u and the (corrected) solution 
U and verify whether the solutions are close to each other in any sense. Such numerical evidence 
could give an idea about the areas of applicability (if any) of the approach. 

The article is organized as follows. In the next section several ways to construct two-scale 
extension for arbitrary initial coefficients aM{x) are presented. The section [3] contains cell problems 
and averaging algorithms from the homogenization theory. Section U consists of numerical results 
in ID for C and extensions (Subsection 14. 2p . and I?2-extension in 2D (Subsection 14. Sp . 

2 Two-scale extensions 

The two-scale extensions (which we numerically investigate in this article) and their properties were 
presented and discussed in [8j . 

The trivial extension is given by a{x,y) = aM{x). More useful extensions can be constructed 
in the following way (assuming that aM{x) is known in a larger domain 17 D $7 in order to avoid 
uncertainties close to 50): 

• we choose e > (small in comparison to the typical size of 17); 
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• for each x G we choose an e-cube Wx with sides aUgned with the coordinate axes, containing 
x: X € Wx- We also assume that Q is large enough: Wx C $7, Vx € {Wx is a cubic 
"Representative Elementary Volume" around x, e is a size of averaging). 

It is reasonable to distinguish two main choices of Wx {C - continuous, P - discrete): 

(C) Wx is an e-cube with the center x; 

(V) Having a partition O = IJ^i i^i H fij = 0, i 7^ j) that each Qj has an e-cube Wj 
(Oj C Wj, Xj is a center of Wj) then for each x G ilj we can define Wx ■= Wj. 

Now we fix X G and construct a(x, •): 

1. d{x,y) = aM{y), y G Wx] 

2. a(x, y) is extended e-periodically in y to the whole space W^; 

3. a(x, y) = a(x, ey) is the two-scale extension. 

Depending on the choice of Wx we have C-extensions and 2?-extensions. 

Remark 2.1. The function a{x,y) is still a two-scale extension if we substitute the item 1. above 
by one of the following more weak requirements: 

• d{x,y) = ajyfiy), y G 0(x) C Wx, where 0(x) is a neighbourhood of x; 

• d{x,y) = aM{y), y = x; 

and let d{x,y) to be free in the rest ofWx- 

This can be used to modify the coefficient near the boundary of Wx e.g. if we want a{x,y) 
to be continuous. The second requirement is so weak that it allows to construct any two-scale 
extension satisfying Def \l.l\ (without saying how to do it). Probably we should have something like 
d{x,y) PS aM{y), y G 0(x) not to loose the relation between OMi') and a(-,-) completely. Anyway, 
we don't consider these possibilities further in this paper. 

The P-extension depends on the choice of {f^j}. If we are going to solve ([2]) using an unstruc- 
tured grid then {^j} could be chosen related to that grid. For example, if we deal with FEM then 
each ilj could be a union of one or more finite elements. Here we will test only one kind of a 
subdivision of into {f^j}, which is more appropriate for solving ([2]) on Cartesian grids: 

Definition 2.1. Let k > 1, e > be given. We divide M'^ into cubes 

□/= (ii/i, (n + l)/i) X ••• X (^irf/i, (id + l)/i), h = e/k, I = {ii, . . . ,id) e Z'^. 

We set fij = n 0,, where is some numeration of those cubes which have a non-empty 

intersection with il, j = 1, . . . , A'x). Wj is a cube with the side e = kh, and the center at the same 
point as the center o/D/q). The D-extension constructed this way let us call a Di^.- extension. 

3 Averaging using two-scale extension 

The sequence ^ is well investigated in the homogenization theory. It is known that the averaged 
coefficient A{x) at x G in the limit problem ([2]) can be calculated via a so-called cell problem. 
Next we remind different formulations of the cell problem applied to the two-scale extensions from 
Section [21 Let us fix an arbitrary x G il. 
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Differential form in W^: 



-V,, • I a 



m 



/ Wj{x,y) dy = 0, Wj{x,y) is 1-periodic in y. 

Y 



(5) 



Differential form in Y. Since we prefer to solve the problems in a bounded domain, we can 
rewrite them in a differential form in a cube Y = (0, 1)°': 



-Vy ■ [a{x,y) {VyWj {x,y) + ej)^ =0 



in Y 



Boundary conditions on Sf,Sj for all i = 1, . . . , d : 

Wj{x,-)\so = Wj{x,-)\sl 

ei ■ {a{VyWj + ej)){x, ■)\so = Ci ■ {a{VyWj + ej)){x, 
fwj{x,y)dy = 



(6) 



where Sf^ = {y & Y : yi = a}. Wj{x, •) is extended periodically in y from Y to 



Variational form: find 'Wj{x, •) G Hp^j.{Y) /R such that 

j^\/y(f>{yfa{x,y)VyWj{x,y)dy = -j^Vy(t>{yfa{x,y)ejdy ycf> € H^,,{Y)/R. (7) 
The averaged coefficient A{x) can be calculated from the solutions wj, j = 1, . . . ,d: 

A.,W=/^.M...)(V,.,(...)+.,)*. (8) 
After solving the solution U could be corrected (see e.g. [I] p.76]): 




Roughly speaking, U approximates u from ([T|) itself, and U approximates the averaged u. 

Due to 1-periodicity of Wj{x,y) in y, we can substitute Y in dS])-® by any other 1-cube 
C = iclnin,ci,in + 1) X • • • X (c^i„,c^i„ + 1). 5f we Can redefine as {y £ C ^ yi = d^^^ + q}. 
For practical purposes it is convenient to take (7 = 1^, where Y^ = {y ^ \ ey G Wx\ for each 
fixed x € O. Then for y G we have a{x,y) = a{x,ey) = aM{Sy). It is also useful for calculating 
the correction ([9]) since x being (always) inside Wx implies x/e € Yx- Thus, we don't need to store 
a(-, •) as a function d x d variables - it is possible to obtain all necessary information directly 
from aM(')- The averaging method is local: the averaged coefficient A{x) and Wj{x, •) depend only 
on the values of aM{-) in Wx, a neighbourhood of x. 

Wx in the C-extension is changing with the point x; the field A{-) is a result of solving the cell 
problems at all points from $7. This is different from the P-extension, where a finite number of 
cell problems has to be solved since Wx is the same in {Wx = Wj). In this case A{x) has a 
constant value in each {J^j}. We note, that the averaged coefficients from both C and P extensions 
coinside at the points xj, the centers of Wj. We also know that ^4(3;) from the C-extension should 
be continuous ([SI Prop. 7.1]). Therefore from a practical point of view these extensions could be 
seen as different interpretations of the coefficient A{x) known at the finite number of points xj: we 
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can treat the data as a continuous or a piecewise constant function. The continuous data can be 
interpolated in space between Xj. The interpolation makes possible the numerical averaging with 
C-extensions. Such averaging needs additional care comparing to the averaging from P-extensions: 
if the distribution of Xj is not dense enough in Q, the interpolated field could be significantly 
different from the exact A{-) (e.g. by missing oscillations). 

Remark 3.1. The Vi~-extension, (0) lead to the averaging algorithm for A{x) proposed in 
fm p. 527] as one of several alternative upscaling procedures using "border regions". The early 
algorithm |^ can be obtained from the T>i-extension (where Wj = Vtj). 

Remark 3.2. The corrected approximation U, calculated from a T>-extension via is not con- 
tinuous. The jumps on J7„ fl o,i^e expected due to the abrupt change of the cell solutions Wj 
when X goes from to 0,m. These jumps are more significant for T)k- extensions with smaller k 
since for large k, Wn and Wm have a large common volume. Here it creates no problem since we 
use only L?' , norms for the comparison U with u. Although if one is interested in fluxes or 
approximations then the correction in the form (0j is probably a bad choice. 

4 Numerical results 

Our main purpose in this section is to solve several model problems ([T|), ([2]) semi-analytically (if 
possible) or numerically and to compare u with U and U . 

4.1 Random number generator 

In most of the numerical examples in this article the coefficient 0^(0 is defined with the help of a 
random sequence {^i}. To generate the sequence {^j} of real numbers we read at each occasion an 
i-th pair (6o, ^i) of bytes from the file [6] (6o> G < 6o> ^ 255) and calculate 



The first five pairs are (34, 178), (52, 184), (220, 178), (237, 13), (19, 247). This approach was chosen 
since it is easy to reproduce the sequence on different computer platforms. 

4.2 ID tests 

In ID we have the following problem 



where a{x) is either the initial coefficient aM{x) or the averaged coefficient ^4(2;). aM{x) is a 
constant in [0,1/4) U (3/4,1] and has oscillations in [1/4,3/4] (see Fig(lH3]). The solution to the 
equation can be written analytically: 



e. = (6o + 6i2«)/(2^^-l) 



eiG[0,l], i = l,2,... 






C 



+ 



F{x) 



(10) 
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where C can be determined from the boundary condition u{l) = u^: 

C = [ —— dx ] { Ur — ui — — -— dx 



aix) J \ Jo a{x) 

Thus, the semi-analytical numerical solutions for initial and averaged problems need only the 
numerical integration. It seems to be more flexible not to consider a(-) in some exact analytical 
form, but to use discretizations of a{x), u{x) on uniform grids. Thanks to one-dimensionality, grids 
with millions of points are available (A^soi ~ number of points). 
The cell problem {x is like a parameter here) 

^ (^a{x, y) (^^^ + l) ) = 0' «^(^' 0) = ^(^' 1) 
also can be solved analytically (up to an additive constant): 

dw{x,y)_ C{x) , ^^,y)=^^,o) + C7(x) 



dy a{x,y) ' ' ' a{x,y) 

where C{x) = ( f} a{x,y)^^ dy^ , since w{x,0) = w{x, 1). The averaged coefficient is 



A{x) = a{x, y) (^^^^lA + l]dy = C{x) 







A{x) 



dy 

dy /^l r dz (l /■™+(^) dz 



lo a,{x,ey)J \e Jq a{x,z)J \£ Jw_{x) aM{z) j 

where Wx = {w-{x),Wj^{x)) = {x{x) — e/2,x{x) + e/2). For the C-extension: x{x) = x. For the 
Pfc-extension: x{x) = h{[x/h\ + 0.5), where [y\ is the largest number from Z: \_y\ < y. 
The correction ([9]) is 

U{x) = U{x) + eU'{x)w{x,x/e), 
where the expression for U'{x) can be found in pop . 

Remark 4.1. The harmonic averaging is used in the finite volume method e.g. for discretizing the 
elliptic operator with discontinuous coefficients fll 



For the ID tests the coefficient a^ix) in 17 = (—1,2) is 

. ^ / 1 xe (-1,1/4)U[XM,2) 

where 

1 (0.1+46^-1) ^ 3 

Xi = -, Xi+i =Xi + e — , XM-l < ^ ^ ^M, 

{^i} is the pseudo-random sequence of numbers, e is either 0.004 (case al, see FigH]), 0.001 (case 
a2, see Figl2]) or 0.00025 (case a3, see Figj3]). The homogeneous boundary conditions g = 
{ui = Ur = 0) are chosen. We use three different r.h.s.: oscillating, constant and discontinuous. 



fix) 



50 sin(30x) case fl 

—4 case f2 

, 4(1(1/2,3/4) - 1(1/4,1/2)) case f3 
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where 'i-(a,b)i^) is a characteristic function of {a,b). 

Let us look at one test more precisely. In Fig iH the averaged coefficient for the case a2, obtained 
from C-extension for e = 0.016, is plotted. From FigJSl where U, U and u are compared, we see that 
the averaging is capable to provide good approximations, and that the correction U approximates 
u with a better quality than U (the later looks more like an average of u smoothing the abrupt 
curve) . 

To estimate quantitatively the quality of the approximation we will use: 



Eo 



\U- 



''"IIl2(o,i)i 



If/- 



" ^IIl°°(o,i)' 



\U- 



■^IIl2(o,i), 



So 



If/- 



■ ^IIl°°{o,i)- 



4.2.1 C-extensions in ID 

In the first series of tests we solve the problems ([I]),© for different aM{-) and /(•) (cases a2fl, a2f2, 
a2f3, alfl, a3fl). The averaged coefficients are calculated from the C-extensions for different e. The 
approximation errors are plotted in Figl6Hl0l In all cases the uniform grids have Ng^i = 8 • 10^, 
16 • 10^ 32 • 10^, 64 ■ 10^ number of points. We can see from the figures, that E2, E^o curves for 
different Nsoi are splitted at the end {e 10~^). Rounding errors and insufficient resolution could 
probably explain this, since the curve obtained on the coarsest grid Nsoi = 8 • 10^ starts to deviate 
first, and the curve from the finest grid N^oi = 64 • 10^ remains longer close to the extrapolated 
line. The numerical results show that smaller e lead to more accurate approximations, and that U 
approximates u better than U does. The curves on some intervals look like straight lines (especially 
E2). The slopes of the lines on the log-log plots give an idea about the order of convergence. 




0.9 



Figure 4: Averaged coefficient A{-) for the case a2 obtained from C-extension for e = 0.016 



4.2.2 C-extensions and D^-extensions in ID 

Calculation of the coefficient A(-) from a C-extension needs high computational resources (due to 
the fine grid), since the fine scale details of the averaged coefficient (see Figjl]) could disappear 
after interpolation of a coarse grid data. Opposite to that, the averaged coefficient from a V- 
extension is free from the interpolation error, and the needed computational resources are limited 
by the particular choice of the extension. Let us compare the qualities of approximation from C, 
Pfc-extensions for k = 1,2,4,8. The grid has Nsoi = 64 • 10^ nodes. From Figllll we see that 
the C-extension provides better U approximations (possibly with higher order of convergence), 
although there is no significant difference when U is concerned. We also observe that the quality of 
approximation from the P^-extensions approach the quality of approximation from the C-extension 
when k increases. 

The (semi)-analytical solutions U, U' , U were used also for the P^-extensions. This means that 
the errors which would appear in practical situation {UhjUj.^ instead of U,U') were excluded here. 
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Figure 5: Comparison of u with U and U for the case a2fl obtained from C-extension for e = 0.016 
4.3 2D tests 

ID case is very favorable for investigations: extremely fine grids and analytical expressions for the 
solutions are available. In 2D we are much more limited in means: we have no analytical solution for 
more or less realistic problem specification, and the finest grid for calculating numerical solutions 
contains only few thousand nodes discretizing OX,OY directions (here the maximum is 4096). 
Appearance of arbitrary directions makes the difference from the ID case. 

A reliable investigation of the C-extension remains practically out of reach here. Thus, we 
restrict ourselves to P^-extensions for k = 2. The extension has one parameter - e. We also use 
the equivalent parameter h = e/k = e/2 emphasizing that the matrix valued coefficient A{-) is a 
piecewise constant function on the /i-grid. A coarser grid cannot resolve the coefficient properly. 

The domain for 2D tests is O = (0, 1)^. The right hand side and the boundary values for ([T]), 
([2]) are fixed for all tests: f{x) = 10 in $7, = on d^l. The coefficients aM{-) are described below. 
We choose only infinitely smooth coefficients to optimize the accuracy of the numerical method on 
available grids. omCO can be naturally extended from Q to any O C M?. 

To solve the 2D elliptic problems with tensor coefficients (fine scale problem ([T]), homogenized 
problem ([2]), cell problems ([7])) we divide the domain (0, 1)*^ by a uniform Cartesian grid into N x N 
squares with the side h = 1/N (/i-grid). All squares are subdivided into two triangles by the same 
diagonal, and the standard finite element method with linear base functions on such triangulation is 
used to solve the problems numerically. The coefficient is forced to have a constant value inside each 
square by taking the value in the center of the square for the whole square (such approximations 
are used for ([1]), ([7]) since the initial coefficients aAf(-) are smooth in our tests). 

The averaged coefficient which is actually used to solve ([2]) numerically is different from the 
exact A(-) due to errors of approximation introduced while solving the cell problems on Nc x Nc 
grids. Let us call it Ah^hd') instead of A{-). The first index h emphasizes that the coefficient is 
piecewise constant on the /i-grid, and the second index he = 1/Nc specifies the discretization step 
used to solve the cell problems. Nc is independent from N = 1/h and should be large enough for 
solving cell problems with enough accuracy. In the tests described below, Nc was usually chosen as 
large as possible under a constrain of reasonable total time of solving A^'^ cell problems on a single 
processor computer. In addition, the grid (NNc/k) x (NNc/k) {k = 2 here) was fine enough for 
resolving all oscillations of aM{-) in In some cases A^^^^ was compared with Ah^2hc^ and the 
solutions of ([2]) with both Ah^hc and Ah^2hc were compared with each other in order to verify how 
the error in A{-) affects the accuracy. 
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10"* 10"' 10"" 10"' 10° 10"* 10"' 10"^ 10"' 10° 10"" 10"' 10"" 10"' 10° 10"* 10"° 10"^ 10"' 10° 



Figure 8: e = 0.001, case £3, C-extension: E2, E^, E2, E^o depending on e 
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10"" 10^ lO"'' 10' 10° 10"' 10"^ 10"^ 10"' 10° 10"* 10"' 10^ 10"' 10° 10"" 10"' lO"'^ 10"' 10° 



Figure 9: e = 0.004, case fl, C-extension: E2, E^o, E2, E^ depending on e 




id" IO"' 10"^ 10' 10° 10" ID"' 10^ id' 10° 10"* 10"' 10^ 10' 1D° ID"" 10 ' 1D^ 10 ' 10° 



Figure 10: e = 0.00025, case fl, C-extension: £^2, E'oo, -^2, E^ depending on e 




10 " 10^ 10^ 10 ' 10° 10 " 10"^ 10^ 10 ' 10° 10 " 10^ 10^ 10 ' 10° 10* 10^ 10^ 10 ' 10° 



Figure 11: e = 0.001, case fl, £^2, -^oo, -E'2, -E'oo depending on e for different extensions: Pi, V2, 
X>4, 'Ds and C 
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The problem ^ with the coefficient we solve numerically on two grids: /i-grid and 

/i/4-grid. The solutions are Uh and Uh^i respectively. Uh is cheap and therefore appropriate for 
solving practical problems, although the (coarsest possible) /i-grid cannot guarantee that Uh is a 
good approximation for U. For example, the difference between Uh and Uh^A is important when 
Ah^hci') has a high contrast. Thus, we need also Uh^i ~ our numerical substitute for U. 

In order to construct the numerical corrections Uh, Uh,A approximating U from ([9]) we need to 
save the solutions of the cell problems. Since the computer memory is also a limited resource, the 
cell problem could be solved on Nc x N^, grid, but saved on Ncs x N^s grid for Ncs < -/Vc- And we need 
to store the values of Wj only at the points which correspond to Vti inside Wi. For example, we can 
choose a priory a set of points {xk} in Vt where we would like to know U , and store the interpolated 
cell solutions from Wi only at the points corresponding to G ilj. The derivatives from [/ in ([9]) 
are approximated in the centers of h squares via central differences and then interpolated in Vt. 
The values in the central differences are either from Uh or from the projection Uh^i to the /i-grid. 

The following relative errors are used to compare the numerical solutions with the reference 
solution: 

E2{y) = \\y - ^^re/||L2(n)/ll^re/llL2(Q), Eoo{y) = \\y - ^^re/ 1| L°° (H) / Ikre/ (f7) i 

where the reference solution Uj-ef is a numerical solution of ([T|) obtained on the finest grid N,i~f,f x 
Nref- Nref is either 2048 or 4096 depending on the intensity of oscillations in aM(')- 

Each Fig ll2|l71120l consists of two subfigures with E2 (left) and -Eqo (right) error functions. On 
each subfigure there are 3 functions: ci{h), C2{h), c^{h). The markers correspond to all test cases. 

ci The curves with square markers represent the functions ci(/i) = E2{uh) for the left subfigure, 
and ci{h) = E^[uh) for the right subfigure, where Uh is the numerical solution of ([I]) obtained 
on the /i-grid without averaging. Uh on the finest grid is the reference solution u^ej and 
therefore the corresponding square markers for E2{uref) = Eoo{uref) = are excluded from 
the curves. 

C2 The curves with circles represent the functions C2{h) = E2{Uh) for the left subfigure, and 
C2{h) = Eoo{Uh) for the right subfigure. 

C3 The curves with point markers represent the functions c^{h) = E2{Uh,A) for the left subfigure, 
and c^{h) = Eoo{Uh,i) for the right subfigure. The averaged coefficient is the same as for C2 
- Ah^hS')i but C3 is different from C2- 



4.3.1 Test with explicitly given coefficient 

In [To] the following coefficient for ([T]) was proposed as a test "without scale separation": 

1 /l.l + sin(27rxi/ei) 1.1 + sin(27rx2/e2) 1.1 + cos(27rxi/e3) 
aM{xi,X2) = - 1^1,1 + sin(2^^2/ei) ^ 1.1 + cos (2^x1/62) ^ 1.1 + sin(27rx2/e3) ^ 

I.l + sin(2vrx2/s4) 1.1 + cos(2vrxi/e5) 2^2^ + ^ 

l.l + cos(27rxi/e4) 1.1 + sin(2^X2/e5) ^ 1 2; ^ 

where £1 = 1/5, £2 = 1/13, £3 = 1/17, £4 = 1/31, £5 = 1/65. 
The curves ci, 02,03 for this test are plotted in FigfT2l 
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Figure 12: ^2(^/^,4), E2{Uh), - Mt, E^{Uh,i), E^{Uh), E^{uh) - right 
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4.3.2 Tests with randomly constructed coefficients 

Let us consider the scalar coefficient aM{x) = 10^'^^^^ where 



S{x) = ^ sin(7ri(xi sin(V'j) + X2 cos(V'i) + 4>i)), ipi = 2ttS,: 



2i-l, <Pi 



26i P 



1=1 



{^j} is the pseudo-random sequence of numbers, the constants m,M 



logio(C) 
M -m' 



iVsin 64 128 256 512 

m -19.7229 -36.1412 -49.6262 -81.8554 
M 22.5351 34.124 51.5507 75.7885 

give approximations to minimum and maximum values of S{x) in 0, respectively. This allows us to 
choose the constant C = 10^ as the contrast for aM(') (C ~ maxx auix)/ min^; aM(a^))- 

We use 4 different coefficients aM(') with different intensities of oscillation: iVgin = 64,128,256,512 
(see Fig J13p . From this series we can observe what happens when aM{') becomes more and more os- 
cillatory, and guess further behaviour towards more realistic situations. One test case {Ngin = 256, 
h = 1/16) is illustrated in Fig ]14ll5] (see also [9], where similar results for another aM(') were pre- 
sented). The curves ci,C2,C3 are plotted in FigfTTl- FigJ^UJ The contrast of the averaged coefficient 
is presented in Fig. [TH 



4.3.3 An interpretation of the presented 2D results for T>2 extensions 

With the help of the information presented in Fig |12l Fig |17fl20l it is possible to estimate the abilities 
of the proposed P2 averaging approach (c2 - practical, C3 - theoretical) in comparison with the 
direct numerical approach (ci). 

For each aM{-) we introduce a level Hres which approximately separates the discretization steps 
{h} into two groups: 1) resolving (/i < Hres) and 2) not resolving (h > Hres) the initial coefficient 
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OAf(-)- Hres is a characteristic value, it is not uniquely defined. We can choose Hres — 1/(2 • 65) 
for the first 2D test (Subsection l4.3.1l) . and Hres ~ 1/^sin for the rest 4 tests (Subsection l4.3.2l) . 

When h < Hj-es^ c\{h) is a monotone increasing (with a constant rate) function of h. In the 
region h > Hres, ci{h) is nearly horizontal since the direct numerical methods fail to approximate 
well problems with rapidly oscillated coefficients until the coefficients are resolved (such behaviour 
is not shown in our figures, except Fig j20|) . 

C2{h), C'i{h) behave in a more complicated way. The upscaling is the most effective for coarse 
grids, h > Hupsc, where C2{h), C3{h) are monotone increasing (with a constant rate) functions of 
h, almost coinsident to each other. To illustrate the choice of Hupsc, we refer to Fig{T2]and Figl20l 
where H^psc — 1/32 and Hupsc — 1/16 respectively. 

When h decreases further, h < Hupsc, the accuracy of the approximation improves but with 
the slowing down rate. The averaging still makes sense, but it is less effective as before. In all 
cases except Figll2( C2 reaches a local minimum at some h = Hacc- Further grid refinement in 
the averaging process gives deterioration in the accuracy. Monotone is a desirable property for 
the 'accuracy vs. discretization size' functions, but unfortunately it is unlikely to hold even for C3 
curve. C2{h) and c^{h) are almost the same for h > H^ev and start to deviate from each other for 
smaller h. This happens since the increasing contrast of A{-) (see Fig |16p prevents the accurate 
solving of ([2]) on the /i-grid. 

We observe that in the region of the resolved aM{-), C2 comes close to ci (with similar slope) 
and possibly crosses it. For small enough h {h = e/2 < Hres) and continuous aui-), the coefficient 
used in cell problems has a small variation. Consequently the averaged coefficient A{-) can be 
seen as a perturbation of aA/(")- Thus, there is no surprise that ([1]), ([2]) after solving on the same 
/i-grid by the same numerical method lead to similar results for h < Hres- Also, we note that it 
is intuitively better to apply a numerical method directly to aMi') than to its perturbation A(-) 
when the grid easily resolves the initial coefficient. This gives some explanation why the averaging 
algorithms rapidly improving at coarse h have to slow down and to 'wait' the direct method. Similar 
behaviour is called "resonance" in the terminology of the multiscale finite element method j^. 

Let us look how the curves change when 0^(0 becomes more and more oscillatory (A'sin increases 
from 64 in Fig(T7]to 512 in Fig |20p : 1) ci moves to the left - Hres decreases; 2) the region where 
the averaging is effective has a tendency to expand - Hupsc, C2{Hupsc) decrease; 3) improving of the 
best accuracy which can be achived on coarse grids (it can be roughly characterized by C2{Hacc) if 
the local minimum exists). 

The quantity 

sup max{^ii(x), A22{x)} 
^ xgn 

inf min{ylii(a;),^22(a:)} 

plotted in Fig[T6]for different e and A'^gin is related to the contrast of A[-). The averaged coefficient 
A{-) is rapidly oscillated when e is small, and A{-) ~ const when e is large. In other words, 
A{xi) ~ A{x2) even if W^^ H = and xi and X2 are far from each other. This could be an 
indication of some statistical properties of our coefficients aM{-), possibly useful for reducing the 
computational cost of the averaging (see the discussion of linear and sub-linear cost of upscaling 
algorithms in [4j.[10j). 

5 Conclusion 

In this article the averaging algorithm for the second order elliptic equation using C and T>k two- 
scale extensions was described in details and applied to several one and two dimensional model 
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problems. Our purpose was to show that there are non-periodic coefficients aA/(') for which the 
standard periodic homogenization together with the two-scale extensions could provide reasonably 
good averaged coefficients. For the test cases we investigated how the quality of the approximation 
depends on the averaging size £, and how the averaged approximations Uh and Uh perform against 
the direct numerical approximation (without averaging) Uh. 

We need to mention that one can construct such initial coefficients aM{-) for which the presented 
here averaging algorithm fails to approximate well on coarse grids. In these cases the averaging 
has no advantage over the direct numerical method. The topic we are planning to address in a 
forthcoming work. 
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Figure 14: Au (left), A12 = A21 (middle), A22 (right) for iVsin = 256, h = 1/16, P2-extension. 
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Figure 15: Comparison of u with U and U on several cross-sections for Nsin = 256, h = 1/16, 
P2-extension. Uref was calculated on 4096^ grid, U on 16^ grid, cell problems on 512^ grids. 
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Figure 16: Contrast of the averaged coefficient depending on h = e/2. The contrasts of qm are 10^ 
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